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Abstract 

Networks have been studied mainly by statistical methods which em- 
phasize their topological structure. Here one collects some mathematical 
tools and results which might be useful to study both the dynamics of 
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agents living on the network and the networks themselves as evolving dy- 
namical systems. They include decomposition of differential dynamics, 
ergodic techniques, estimates of invariant measures, construction of non- 
deterministic automata, logical approaches, etc. A few network examples 
are discussed as an application of the dynamical tools. 

1 Introduction 

When modeling extended complex systems, the network concept appears quite 
often. The metabolic processes of living beings are a network with the substrates 
as nodes, linked together whenever they participate in the same biochemical 
reaction. Protein-protein as well as gene expression and regulation are also 
networks. Social, economic and political networks are the backbone of human 
society, the internet is a network, etc. [Albert & Barabasi, 2002] [Dorogovtsev 
& Mendes, 2003]. Most studies deal with networks as statistical objects, with 
extensive use of the tools of statistical mechanics [Pastor-Satorras et al., 2003]. 
Much less attention has been paid to the dynamical phenomena taking place in 
the networks or to the behavior of the evolving networks as dynamical systems. 

For several decades an intensive effort has been dedicated to the study of 
low-dimensional dynamical systems, leading to an extensive body of rigorous 
results [Katok & Hassclblatt, 1995]. This exploration is still proceedings at a 
good pace with exciting new results, for example, in the dimension theory of 
dynamical systems [Pcsin, 1997] [Barreira, 2002] and non-uniform hypcrbolicity 
[Bonatti et al., 2003]. However, the main challenge for physical applications lies 
on extended systems, in particular on the understanding of the dynamics leading 
from microscopic laws to global patterns [Cross & Hohenberg, 1993]. A large 
amount of numerical work has been done on the dynamics of these systems which 
led to the identification and classification of typical patterns, spatio-temporal 
chaos, statistical properties and multistability [Kaneko, 1993] [Kaneko & Tsuda, 
2000] [Boccaletti et al, 2001] [Boldrighini et al, 2000]. Rigorous results are few, 
except for regular coupled map lattices [Bunimovich & Sinai, 1988] [Coutinho 
& Fernandez, 1997a, 1997b] [Jiang & Pesin, 1998] [Afraimovich & Fernandez, 
2000] [Gielis & MacKay, 2000] [Fernandez & Guiraud, 2004]. 

For non-uniform coupling structures there are much less results, statistical 
mechanics tools being used mostly to characterize the topological features of 
static and evolving networks. Of course the topological structure of the network 
is very important for the dynamics. Form affects function and topology controls 
the rate at which information or diseases propagate [Boots & Sasaki, 1999] 
[Keeling, 1999] [Pastor-Satorras & Vespignani, 2002] [Lloyd & May, 2001], the 
robustness under attack and failure [Albert et al., 2000] or the adaptation and 
learning processes that take place [Araujo & Vilela Mendes, 2000]. 

The main purpose of this review is to provide a toolkit for the treatment 
of networks (both regular and irregular) as dynamical systems. Results from 
differential dynamics and ergodic theory will be presented. To deal with a 
network as a dynamical system, three main problems have to be addressed. 
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First, how to characterize the dynamical behavior of the "agents" sitting on the 
nodes with interactions defined by the link structure of the network. Second, if 
the network topology is not fixed either because the links are changing in time 
or because the network itself is growing, to characterize the network evolution as 
a dynamical system. Finally to characterize the interplay between the network 
topology and the nature of the dynamics. 

Section 2 concerns the description of network dynamics by global functions, 
limit cycles and multistability, with some examples illustrating both dynamics 
on a network and network evolution as a dynamical system. Section 3 discusses 
ergodic techniques (old and new) and their role in networks. Section 4 discusses 
the logic approach to network dynamics. On the important problem of topology 
versus dynamics, I will touch only briefly in the ergodic section, referring for a 
particular detailed example to [Araujo et al, 2003]. 

In almost all cases the reader should refer to the original papers for proofs 
and further developments. Occasionally, as for example in the family of ergodic 
parameters or in the relation between synchronization and Lyapunov spectrum, 
somewhat more detail is given to the exposition. This does not mean in any 
way that I consider this to be a more important topic than the other inter- 
esting results of many authors that are mentioned. If more detail is included 
it is because I believe it to be a new result not included before in any other 
publication. 

2 Differential dynamics tools 

2.1 Describing dynamics by global functions 

In many networks, the node dynamics may be modelled by ordinary differential 
equations of the form 



For a neural network, the x\s might be firing rates and the W(jS synaptic in- 
tensities [Grossberg, 1988], for a genetic regulatory system [Tyson & Ohtmer, 
1978] [de Jong, 2002] the variables Xi would code for the concentrations of RNA, 
proteins or other metabolic components and Wij for the production constants 
(measuring the strength of j on i), f (•) being the regulation function and —"faXi 
a degradation term, etc. 

2.1.1 Symmetric systems 

Eq.QJ is a particular case of the Cohen-Grossberg form [Cohen & Grossberg, 
1983], used by these authors to describe continuous-time neural networks, 



dxi 
~dt 



Xi (x) = OLi + ^ W ijf ( X j) - li X i 



(1) 





3 



Cohen and Grossberg proved that, for the symmetric case (Wij = Wji), the 
following function 

n ~ Xi -i n 

y (*o = - E / fci(6)/i(6)«ie< + 2 E (3) 
»=i j.fe=i 

is a Lyapunov function, that is 

±V&)<0 (4) 

along the orbits if ai(xi)f i (xi) > 0- Hopfield's [Hopfield, 1984] "energy" func- 
tion is a particular case of this result. 

The existence of a Lyapunov function is a useful device to characterize the 
asymptotically stable states of the network or for the synthesis of networks with 
a desired number of stable asymptotic solutions [Cohen, 1992]. 

In the case of symmetric connections the continuous-time result of Cohen 
and Grossberg has been extended to a class of discrete-time systems ([Fogelman 
Soulic et at, 1989] and references therein). For non-symmetric connections of 
particular form, namely 

HjWij = ii, W (5) 
/ii > 0, and time evolution of the connection strengths of Hebbian type 

Wij = -7y Wij + fi (xi ) fj (xj ) (6) 

in Schiirmann [1989] or WijWji > and J\ c Wij = Y\ c Wji along every cycle 
in Fiedler & Gedeon [1998], Lyapunov functions may also be constructed. 

2.1.2 General systems 

The Cohen-Grossberg result has been generalized for arbitrary w'^s in Vilela 
Mendes & Duarte [1992], namely given 

= wjp + w}? 
w^p = \{Wa + Wn) 

wi A) = I (W {j -Wj^ (7) 

one has the following 

Theorem [Vilela Mendes & Duarte, 1992] If a^Xi)/^ (xj) > Vx,z and the 
matrix has an inverse, the vector field *i in Eq.© decomposes into one 

TT ■ • • (G) • (H) 

gradient and one Hamiltoman component, Xi—Xi + x\ , with 

f'.( Xl ) dxi ~ l^jy*A x ) dx 3 



X 



{H) = v „,„,;,,; (x) aj (xj)^ v r„ ( .r^ 
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and 

WiiW = -a i (x i )- 1 {W( A )-i)..(x)a j (x j )-i 

{^Ylij^ijOJjk — Sikj- 9ij( x ) and u>jk(x) are the components of the Riemannian 
metric and the symplectic form. 



The conditions on cii(xi), fi(xi) and w\f insure that g is a well defined 
metric and that u> is non-degenerate. 

The decomposition ijSJ) is useful, for example, on the design of oscillatory 
networks and on the study of gated learning rules [Howse et at, 1996]. The 
nature of the dynamics in the network will depend on the relative strength of 
the gradient and the Hamiltonian components. Howse et al. [1996] propose to 
measure this relative strength by comparing dv J t - with ^S-. However, these 
quantities vary in space and time and it is the compensation of the two effects 
that in particular regions of phase space lead to the attractors of the dynamics, 
for example to limit cycles (see below). 

The identification, in the differential system of just one gradient and 
one Hamiltonian component, with explicitly known potential and Hamiltonian 
functions, is a considerable simplification as compared to a generic dynamical 
system. For a general dynamical system a representation by one or two functions 
is possible only locally [Vilela Mendes & Duarte, 1981] and explicit forms for 
the functions are not easy to obtain [Abarbanel & Rouhi, 1987] [Crehan, 1994]. 
Global decomposition for general dynamical systems require one gradient and 
n — 1 Hamiltonian components [Vilela Mendes & Duarte, 1981], namely 

- £ w w&r+EE W (*>) l3 ( 10 ) 

3 3 fc=lj=l 3 

{ (x) } being a set of canonical symplectic forms adapted to each Hamiltonian 
component. This result is a generalization to n dimensions of the 2-dimcnsional 
result of Roels [1974]. The first term in (|10fl is the dissipative component and 
the second one corresponds to a volume-preserving dynamical system. 

The above results lead to a convenient characterization of dynamical systems 
of type or @ ■ For the symmetric case the existence of a Lyapunov function 
guarantees global asymptotic stability of the dynamics. However not all vector 
fields with a Lyapunov function are differentially equivalent to a gradient field. 
Therefore the fact that a gradient vector is actually obtained gives additional 
information, namely about structural stability of the model. A necessary con- 
dition for structural stability of the gradient vector field is the non-degeneracy 

,fV <-■>"> 



of the critical points of V 1 -^, namely det 

dxi 



dxidxj 



7^ at the points where 

— t — j 

— 0. In a gradient flow all orbits approach the critical points as t — * oo. 
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If the critical points are non-degenerate, the gradient flow satisfies the condi- 
tions defining a Morse-Smale field, except perhaps the transversality conditions 
for stable and unstable manifolds of the critical points. However because Morse- 
Smale fields are open and dense in the set of gradient vector fields, any gradient 
flow with non-degenerate critical points may always be (^-approximated by a 
(structurally stable) Morse-Smale gradient field. Therefore given a symmetric 
model of the type (0), the identification of its gradient nature provides a easy 
way to check its robustness as a physical model. 

Although Lyapunov functions may in some cases be constructed for discrete- 
time systems [Fogelman Soulie et a/., 1989], the natural functional representation 
of maps is through generating functions. This is well known for canonical maps 
of symplectic manifolds [Amiet & Huguenin, 1980] and has been generalized in 
[Vilcla Mendes, 1986] for noncanonical maps. 

The representation of network dynamics by global function applies to neural 
networks of several types [Grossberg, 1988], to more general networks [Chua, 
1988a, 1988b] and, in view of an established correspondence [Doyne Farmer, 
1990], to a large range of connectionistic systems. 



2.2 Cycles 

Existence of limit cycle oscillations in networks is an important issue [Gouze, 
1998] [Plahte et ai, 1995] [Snoussi, 1998]. The decomposition theorems provide 
a tool to look for candidate orbits with limit cycle properties. Many years ago 
Pontryagin [1934], studying small perturbations of Hamiltonian fields on the 
plane 

• 9H . , . • dH 

x=—+eA(x,y,s), V= -— + sB (x,y,e) (11) 

introduced the notion of generating cycle 7 (c) , lying on a level curve H = c, 
when the perturbed equation has a cycle that depends continuously on e, for 
small |e|, and tends to 7 (c) when e — ► 0. Pontryagin's result states that if 7 (c) 
is a generating cycle, then 



I (c) = / (Bdx - Ady) = (12) 

Jj(c) 

the integration being along 7 (c) at e = 0. 

Further results on the existence of cycles were later proved both for weakly 
coupled oscillators and for more general systems with parametrized families 
of solutions (see [Hoppensteadt & Izhikevich, 1997], chapter 9 and references 
therein). A generalization of Pontryagin's result to dynamical systems with 
constants of motion [Duarte & Vilela Mendes, 1983], leads to a necessary con- 
dition for the existence of a cycle using the decomposition in (|1C)|> , namely 



( vh w • w) + w (*o ( vh(1) ■ VH(k] ) = ( 13 ) 
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the integration being along a closed level curve 7$ of Hi. 

A similar result holds for discrete-time maps which belong to a differentiablc 
arc with constants of motion [Vilela Mendes & Duarte, 1982]. A constant of 
motion for a map / defined on a manifold M is a differentiable function $ : 
M — > R such that for some orbit 7, $ o 7 ^constant. It generalizes the notion 
of first integral which would require this to hold for all orbits. A family of maps 
f £ is called a differentiable arc with constants of motion if (i) each f £ has a 
constant of motion $ e for some orbit j e ; (ii) The constant of motion $0 of /o is 
a first integral in a neighborhood of 70 ; (iii) the maps e — > f e , e — > j e , e — ► $ e 
arc differentiable. Then 



iVo being the period of the orbit 70 . 

Both and (|14(l give only necessary equations for the existence of limit 
cycles in the composite dynamics. Nevertheless they are useful tools to identify 
limit cycle candidates. Sufficient conditions may also be obtained in particular 
low-dimensional cases [Vilela Mendes, 1988, 2000a]. 

In the same way as the Hamiltonian components of the dynamics provide 
a tool to look for limit cycles, the stationary points of the gradient potential 
provide information on the multistability of the dynamics and the nature of 
their basins of attraction. It is also a tool for the construction of the invariant 
measures of the dynamics (see below). 



Existence of multiple stable states with distinct basins of attraction plays a sig- 
nificant role in the dynamics of networks, for example in those associated to 
the basic processes of life. A genetic regulatory network with different stable 
patterns of gene activation explains the emergence of different phenotypic ex- 
pressions in the absence of genetic differences [Laurent & Kellershohn, 1999]. 
Examples are also found in population dynamics [Henson, 2000] , neural dynam- 
ics [Skarda & Freeman, 1987] [Freeman, 1992], geophysics, etc. [Vilela Mendes, 
2000a]. Extensive numerical work has been done on classifying different pat- 
terns of multistability and their basins of attraction (see for example [Wuensche, 
2002]). Here, I will concentrate on the dynamical mechanisms leading to the 
existence of multiple attractors. Although some of these results are only rigor- 
ously proven for low dimensional systems, their relevance for high dimensional 
systems is to be expected. 

(i) In the particular case of networks with symmetric connections the attract- 
ing critical points of the potential function are stable asymptotic states of 
the dynamics. 

(ii) Homoclinic tangencies (the Newhouse phenomenon) 

Contrary to earlier conjectures that generic systems might only have finitely 
many attractors, Newhouse [1970, 1974, 1979] proved that a class of diffeomor- 
phisms in a two-dimensional manifold has infinitely many attracting periodic 



N -l 




(14) 



2.3 Multistability 
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orbits (sinks), a result that was later extended to higher dimensions [Palis & 
Viana, 1994]. For two-dimensional manifolds the result is: 

Theorem (Newhouse, Robinson [1983]) Let / M be a C 3 map in a 2-dimensional 
manifold with C 1 dependence on \i and |det(T a /™ )| < 1 and let the non- 
degenerate homoclinic tangency be crossed at non-zero speed at \i = [io- Then 
for Ve > 0, 3(^i,yU 2 ) C (/xo 5 Mo + £) and a residual subset J C (^1,^2) such that 
for n E J , / M has infinitely many sinks. 

Models of such diffeomorphisms were constructed by Gambaudo and Tresser 
[1983] and Wang [1990] proved that the Newhouse set has positive Hausdorff 
dimension. After these results, intense research followed on the unfolding of 
homoclinic tangencies and an essential question was whether, in addition to 
infinitely many sinks, there would also be infinitely many strange attractors 
near the homoclinic tangencies. The question was positively answered by Colli 
[1998]. The main result is: 

Theorem (Colli) Let /o E Diff°°(M) be such that /o has a homoclinic 
tangency between the stable and unstable manifolds of a dissipative hyperbolic 
saddle p n . Then, there is an open set f2 C Dif f°°{M) such that 

(a) /„ G H 

(b) there is a dense subset D C fi such that for all / E D , / exhibits 
infinitely many coexisting Hcnon-like strange attractors. 

Having established the existence of infinitely many sinks and infinitely many 
strange attractors near homoclinic tangencies, a question of practical impor- 
tance is the stability of the phenomenon under small random perturbations of 
the deterministic dynamics. It turns out that the answer to this question is neg- 
ative. Therefore under small random perturbations only finitely many physical 
measures will remain. 

Theorem (Araujo [2000]) Let / : M — > M be a diffcomorphism of class 
C r ,r> 1, of a compact connected boundaryless manifold M of finite dimension. 
If / = fa is a member of a parametric family under parametric noise of level 
e > 0, that satisfies the hypothesis: 

There are K E N and £ > such that, for all k > K and x E M 

(A) f k (x,A) DB fc (xU ) ; 

(B) f k {x,v°°) « m ; 

then there is a finite number of probability measures fii, ■ ■ ■ fii in M with 
the properties 

1. fii, ■ ■ ■ m are physical absolutely continuous probability measures; 

2. supp/Uiflsupp/ij for all 1 < i < j < I ; 

3. for all x E M there are open sets V\ = V\(x), ■ ■ ■ ,Vi — Vi(x) C A such 
that 

(a) Vi n Vj = 0, 1 < i < j < I ; 

(b) (A\ (Vi U • • • U Vi)) = ; 

(c) for all 1 < i < I and v 00 — a.e. t E Vi we have 

hm - V^ffM)) = / 
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for every <j> G C(M,R). Moreover the sets V\(x),--- , Vi(x) depend continu- 
ously on x e M with respect to the distance d v (A, B) = u 00 (A/SB) between 
mod subsets of A. 
(hi) Small dissipative perturbations of conservative systems 
Conservative systems have a large number of coexisting invariant sets, namely 
periodic orbits, invariant tori and cantori. By adding a small amount of dissipa- 
tion to a conservative system one finds that some of the invariant sets become 
attractors. Of course, not all invariant sets of the conservative system will sur- 
vive when the dissipation is added. However, for sufficiently small dissipation 
many attractors (mainly periodic orbits) have been observed in typical systems. 
Poon & Grebogi [1995], Feudel et al. [1996] and Fcudel & Grcbogi [1997] have 
extensively studied these effects in the single and double rotor, the Henon map 
and the optical cavity map. They find a large number of attractors for a small 
amount of dissipation, in particular in the double rotor map. The large number 
of coexisting stable periodic orbits has a complex interwoven basin of attrac- 
tion structure, with the basin boundaries permeating most of the state space. 
The chaotic component of the dynamics is in the chaotic saddles embedded in 
the basin boundary. The systems are also found to be highly sensitive to small 
amounts of noise. The problem of migration between attractors and their stabil- 
ity in multiple-attractor systems has also been studied by other authors [Weigel 
& Atlee Jackson, 1998] [Kaneko, 1997] [Dutta et at, 1999]. 

Rigorous results may be obtained in particular cases using the ideas of de- 
formation stability [Duarte & Vilela Mendes, 1983] [Vilela Mendes & Duarte, 
1982] [Vilela Mendes, 1988] [Lima & Vilela Mendes, 1988]. For example, let an 
e— family of maps be 

x' = bx + y + f(x, y, e) ^ 

y = y + g(x,y,s) 

with f(x, y, 0) = g(x, y, 0) = For e = the map has marginally stable periodic 
orbits of all periods. Under perturbation some of the orbits become stable ones. 

Theorem [Vilela Mendes, 2000a] If / and g are jointly C 2 in (x,y,e) with 
f(x, y, 0) = g(x, y, 0) = 0, there is an e such that for |e| < \e\ an interior orbit of 
period p of the unperturbed map becomes a stable orbit of the perturbed map 
if and only if: 

(1) P -£d e g(x ( n\y^\0)\ £ =o=0 

n=0 

(2) ed £ E \d x g(x { °\y(°\ e) + (1 - b)d v g(xg\yW,e)\ \ E=0 < 

n=0 L > 

(iv) Coupled dynamical systems near period-doubling accumula- 
tion points 

An example is a system of two coupled quadratic maps [Carvalho et al., 
2001], 

xi(t + 1) = 1 - /U* ((1 - e)x 1 (t) + ex 2 {t) f , . 

x 2 {t + l) = l- ^ {ex l (t) + {l-e)x 2 {t) f 1 ' 
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with x € [—1,1], and /x* = 1.401155... , which is the parameter value of the 
period doubling accumulation point. The result is that for any N there is a 
sufficiently small e (N) such that there are N distinct stable periodic orbits. 



2.4 Network examples 

2.4.1 An excitatory-inhibitory network 

Excitatory-inhibitory networks exhibit a rich variety of activity patterns. They 
have been identified as underlying several biological processes like image segmen- 
tation, sleep rhythms, control of movement in the basal ganglia, etc. [Terman, 
2002]. Here one considers a simple network with two populations, one composed 
of excitatory and the other of inhibitory cells. The equations of motion are 



E-f,v + i /(*;)+ A? (*) 



i = l,. 

i = N - 



,2N 



(17) 



The first population (i = 1, • • • , TV) receives a time-dependent driving external 
signal g (t) and inhibitory inputs from the second population. The second pop- 
ulation receives excitatory inputs from the first population and has a decay rate 
—7. The activation function is 



/ ( x ) = 1 - exp {-fix) 



(18) 



In terms of the global functions described before, the dynamics has a simple 
form 

• _ dV ^ dH 

dxi ^ ,J 



dxj 



(19) 



with 



r being the matrix 1 



V 
H 



( 




1 
1 

V 1 



= -&K*)I>i + ?I>? 

= '>v,r f^xK 



-1 -1 
-1 -1 



-1 -1 

1 
1 



1 







-1 \ 

-1 

-1 





j 



(20) 



(21) 



Except for the contribution of the driving signal and the decay constant, 



the dynamics is Hamiltonian, the symmetric connection coefficients 



w 



X s ) 







vanishing identically. Therefore, the reaction of the network to an external signal 

x the degeneracy of the symplectic form is lifted by a particular choice of coordinates 
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100 200 300 400 500 600 

Figure 1: Dynamics of the excitatory-inhibitory network (|17|l (a = 1,(3 = 
0.1,7 = 0-1) f° r a square- wave signal g (t) 

is simply a damped oscillation. This is illustrated in Fig.l where, starting from 
a random initial condition and after a relaxation period, a square wave g (t) with 
zero baseline is applied to the network. The base level of each unit depends on 
the initial conditions. 

Notice however that, if the agents in the network are biological cells, their 
activation cannot be negative. Then, the dynamical system should be modified 

from Xi= Fi (x) (where Fi (x) is the right hand side of H17|0 to 

*i= Fi (x) ■ OR (sign (xi) , sign (Fi)) (22) 
OR being the logical function 



+ - 



1 


1 


\ — '■ 






(23) 



In this case, because the right hand side is no longer of Cohen-Grossberg 
form, the decomposition l|19l) — l|2U|l gives only qualitative information on the 
dynamical behavior of the system. This is illustrated in Fig. 2 where the same 
square- wave driving force, as in Fig.l, is applied to the modified network. 

2.4.2 A gene regulation network 

The p53 gene was one of the first tumour-suppressor genes to be identified, its 
protein acting as an inhibitor of uncontrolled cell growth. The p53 protein has 
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0.2 r- 

0.18 - 
0.16 - 
0.14 - 



0.12 - 




100 200 300 400 500 600 



Figure 2: Dynamics of the excitatory-inhibitory network l|22ll (a = 1,(3 = 
0.1,7 = 0-1) f° r a square- wave signal g (t) 

been found not to be acting properly in most human cancers, due either to 
mutations in the gene or inactivation by viral proteins or inhibiting interactions 
with other cell products. Although apparently not required for normal growth 
and development, p53 is critical in the prevention of tumour development, con- 
tributing to DNA repair, inhibiting angiogenesis and growth of abnormal or 
stressed cells [May & May 1999] [Vogelstein et ai, 2000] [Woods & Vousden, 
2001] [Taylor et al, 1999] [Vousden, 2000]. In addition to its beneficial an- 
ticancer activities it may also have some detrimental effects in human aging 
[Sharpless & DePinho, 2002]. 

The p53 gene does not act by itself, but through a very complex network 
of interactions [Kohn, 1999]. Here I will discuss a simplified network, which 
although not being accurate in biological detail, tends to capture the essential 
features of the p53 network as it is known today. In particular, several different 
products and biological mechanisms are lumped together into a single node 
when their function is identical. The network is depicted in Fig. 3. The arrows 
and signs denote the excitatory or inhibitory action of each node on the others 
and the letters b, g,c,r,p,m,a denote their intensities (or concentrations). 

The p53 protein is assumed to be produced at a fixed rate (k p ) and to be 
degraded after ubiquitin labelling. The MDM2 protein being one of the main 
enzymes involved in ubiquitin labelling, the inhibitory node (to) is denoted 
MDM2. There is a positive feedback loop from p53 to MDM2, because the 
p53 protein, binding to the regulatory region of the MDM2 gene, stimulates the 
transcription of this gene into mRNA. 
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Blood vessel 
formation (b) 



Oncogenes (g) 



Cell growth (c) 



Abnormality 
detection (a) 



pl4 ARF (r) 



p53 (p) 



MDM2 (m) 



Figure 3: A simplified p53 network model 



Under normal circumstances the network is "off" or operates at a low level. 
The main activation pathways are the detection of cell anomalies (a), like DNA 
damage, or aberrant growth signals, such as those resulting from the expression 
of several oncogenes (the pl4 ARF pathway, r) . They inhibit the degradation of 
the p53 protein, which may then reach a high level. There are several distinct 
activation pathways. For example, in some cases phosphorylation of the p53 
protein blocks its interaction with MDM2 and in others it is a protein that binds 
to MDM2 and inhibits its action. However, the end result being a decrease in 
the MDM2 efficiency, they may both be described as inhibitory inputs to the 
MDM2 node. 

The p53 protein controls cell growth and proliferation, either by blocking 
the cell division cycle, or activating apoptosis or inhibiting the blood-vessel 
formation (b) that is stimulated by several tumors. In our simplified p53 network 
all these effects are coded on the following set of equations 

~£ = kp Wp m f m (to) 

% = W mp f p (p) - W mr f r (r) - W ma a - i m m 

§ = W bc fc (c) - W bp fp (p) (24) 

I = W cg g + W cb f b (b) - W cp fp (p) 

It = W rg g — 7r r 
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One should note that an increased level of cellular p53 is not by itself sufficient 
for it to become a transcriptional activator controlling cell growth. Conforma- 
tional changes of the protein are also needed which are stimulated by the acti- 
vation pathways or may be therapeutically induced. Also some viruses produce 
proteins that inactivate p53. All this means that in reality some of the coupling 
constants in Eas.(|24p. (for example W cp ) may also be dynamical variables. 




-0.1 1 ' ' ' ' 1 0.1 1 ' ' ' ' 1 

02468 10 02468 10 



Figure 4: Time evolution of the network (|28|l 

The regulation functions / (•) arc positive nonlinear functions with a thresh- 
old and a saturation level. By shifting variables to compensate for thresholds 
and rescaling the coupling constants they may be normalized by the coefficient 
of the linear part, that is 

fi (as*) = Xi + ■ ■ ■ (25) 

With a rescaling of p, m, b, c, r and redefinition of the constants we may 
consider 

k p = W mp = W bc = W cg = W rg = l (26) 
Furthermore, from the last equation in (|24(l 

r(t) = -g+(r(0)-^-)e-^ r (27) 

Replacing r by its steady state value g/~f r , and rescaling W mr we are left with 
^ = 1 - W pm f m (m) 

% = f P (p) ~ W mr g - W ma a - j m m . , 

| = fc(c)-W bp f p (p) { *> 

1 = g + w cb f b (b) - w cp f p (p) 
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Figure 5: Same as Fig. 4 with different initial conditions 

a system of four dynamical variables and two control parameters g and a. 
Using the dynamics decomposition discussed in Section 2 one obtains 

V {s) = -J P f' P (0dC-J C f'c(0^ + J m (W mr g + W ma a + lm 0f m (0di 

+|E 1 , =P , m , ii c< ) /.(^/ J fe) 

H = J2 Xt=p , m .b.c I ' fi (0 d Z 

(29) 

9a = tttt^ 1 v = ( 3 °) 









pm 


i (WW - 1) 


k (W pm + 1) 


pb 






pc 


l,w cp 




mb 








mc 








be 


-|(W ci , + l) 


£ (^ cb - 1) 



The reaction of m and c to external stimuli (g and a) and the production rate 
of p are coded on the first three terms of the potential function V^ s \ For 
coupling constants of order unit, one sees from (|31|l the existence of a damped 
Hamiltonian oscillation for the p—m system, and a dangerous runaway behavior 
of b — c arising from its dominantly gradient dynamics. The action of p on b and 
c is of mixed gradient-Hamiltonian type. Hence, from inspection of the nature 
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of the global functions describing the dynamics, one concludes that (at least 
in this model) the controlling action of p53 may only be effective in particular 
circumstances. That is, it will depend on the initial conditions. This conclusion 
is now checked by a detailed study of the solutions. 

Consider first the linear approximation to the system. The solutions are, for 
the p — m system 

p(t) = P+p'{t) (32) 
m(t) — m +m' (t) 

with 

P 
m 



W mr g + W ma a + $ 



(33) 



p> (t) = e-^l 2 | [p (0) - p\ cos at + I hf- (p (0) - p) - W pm (m (0) 

TO ' (t) = e-^/ 2 I (m (0) - ^) cosat + ± (p (0) - P -?f (to (0) - ) 

(34) 

and a = ^/Wpm - 7 2 /4- 

As expected, one sees a damped oscillatory behavior of the p — m system 
and, in the absence of stimuli (a = g = 0) the p level is small and controlled by 
the degradation of to. 

For b and c one now obtains 



Ht)\_ W bp P _ ,* f Wbpp > (T) n / c (0) - W 6p p 



sin ai 



sin at , 



where A is the matrix 



(35) 



This matrix has eigenvalues ±y/W c b implying that b(t) and c(t) are going 
to have terms proportional to exp {ty/W c b) and exp {—ty/W c b) ■ Hence p (p53) 
will only have a controlling effect on cell proliferation if the coefficient of the 
exponentially growing terms becomes negative. Multiplying (|35(l on the left by 
>W7h 1 



the matrix \ [ , T c Tr „ I that diagonalizes A one obtains the coefficient of 
2 V -VWcb J 

the exponentially growing term 

(37) 



B{t) = ^ e(0)-IT^ +j l(0)-J^2 



Jo e^-)^ (^F^ bp + i^cp) P' (r) dr 
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Figure 6: Typical equilibrium configuration of network connections evolved ac- 
cording to Eqs. (a = 1, = 0) 

The conclusion is that control of cell proliferation is obtained only if 3t 
such that B (t) < 0. Therefore it depends strongly on the initial conditions. 
This conclusion, inferred both from the dynamical decomposition and the linear 
approximation is borne out by simulation of the nonlinear problem. Figs. 4 
and 5 show two time evolutions of the equations l|24|l with / (x) — tanh (x) , 
Wpm = W mr = W ma = W bp = W cb = W cp = 1, j m = 0.01, g = a = 1 and 
the vector field Xi= Fi truncated to Xi= Fi-OH(sign(xi) , sign(Fi)), because 
concentrations cannot become negative. The behavior depends strongly on the 
value of the initial conditions. In conclusion, the implication (of the model) is 
that unless p53 starts acting soon enough its action is useless and other means 
have to be used to control cell proliferation. 

2.5 Evolving networks 

In many networks found in Nature, as important as the structure of the net- 
work, is the path that the network took to reach that final state. Social or eco- 
nomic networks, industrial, transportation and communication networks, eco- 
logical webs, biological networks, all are examples of evolving networks. In 
many cases their complex structure is a simple consequence of the principles of 
their growth. Several network growth schemes have been studied (see Albert & 
Barabasi [2002] , Dorogovtsev & Mendes [2003] and Pastor-Satorras et al. [2003] 
for reviews). 

Network evolution occurs either by the addition or elimination of interactions 
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Figure 7: Degree distribution of the network in Fig. 6 



between existing nodes or by the addition of new nodes. In both cases, network 
evolution may be looked at as a dynamical system in the space of network 
connections. In the case of growing networks, this dynamical point of view 
may also be used by considering the evolution from zero of previously vanishing 
connections. 

This dynamical approach will be explored here. Using the global function 
description, discussed in Section 2.1, two types of evolving networks will be 
considered. The simplest situation occurs when the dynamics of the connections 
is derived from a potential. In this case, exact expressions for mean values and 
invariant measures may be obtained. 

Consider 



V 1 ({W}) = a'£w? j (W ij 

i<j 



with the network evolving according to 



dt 



dV x 
dW~ 



(38) 



(39) 



When d/0 and /3 = 0, the connections evolve either to zero or to one, de- 
pending on the initial conditions. Therefore the network (with TV nodes), as 
a dynamical system, is a multistable system with 2 N ( N ~ 1 ^ 2 different equilib- 
rium points. A typical configuration, obtained from random initial conditions, 
is shown in Fig. 6 (N = 100) to which corresponds the degree distribution shown 
in Fig.7. 
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Figure 8: Same as in Fig. 6 with a = 1 and (3 = 0.003 



When [3 =/= Q the behavior is quite different, as shown in the typical config- 
uration of Fig. 8 and degree distribution Fig. 9. The degree Ki of a node i is 
defined to be 



holding for all intermediate values of Wij . 

One sees that for (3^0 some nodes are more connected than others. 
V\ ({W}) with (3 7^ is a model for preferential attachment. 

It is not be practical to obtain mean values and distributions directly from 
simulations. This being a multistable system many different simulations with 
well distributed initial conditions would be required to obtain accurate values. 
However, in this case, exact expressions may be obtained from the unique in- 
variant measure for the system with small random perturbations, as discussed 
in Section 3.3 



V 2 ({W}) ^a^Wl 3 (W v - 1) 2 +(3J2 E T—\ ( W ik + W lk) ((Wik if + (W jk l) 2 ) 

(42) 



For [3 a typical configuration is shown in Fig. 10. The main feature is the 
correlation between node connections. For a = 1, (3 = 0.05 and N = 100 
the sum of correlations between the node connections is around 20 whereas for 
a = 1, (3 = it is w 4.5. In conclusion, V% {{W}) is a model for (approximate) 
node replication. 




(40) 




(41) 



As a second example consider 
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Figure 9: Degree distribution of the network in Fig. 8 

3 Ergodic tools 

Topological and differential notions provide useful characterizations of the over- 
all structure of phase space. However, what is more important for the applica- 
tions is the dynamics in the phase space regions most frequently visited by the 
system. This is provided by the ergodic theory, in particular by the classification 
of invariant measures and their characterization by ergodic parameters. 

Let a dynamical system evolve on the support of a measure /i which is left 
invariant by the dynamics. An ergodic parameter If (/u), characterizing the 
measure, is obtained whenever the following limit 

1 T 

MM)= lim = (/»*„) (43) 

n=l 

exists for ^—almost every x - For continuous-time dynamics / denotes the 
time-one map. 

3.1 Lyapunov and conditional exponents 

Lyapunov exponents are the most widely used ergodic parameters. More re- 
cently conditional exponents have also been proposed as an useful characteriza- 
tion of the dynamics. 

Let / : M — > M , with M C R m , [i a measure invariant under / and £ 
a splitting of M induced by R k x R m ~ k . The conditional exponents are the 
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Figure 10: Typical equilibrium configuration of network connections evolved by 
the potential V 2 ({W}) (a = 1,(3 = 0.05) 

eigenvalues and ^ m of the limits 

lim (D k f n *(x)D k f n (x))^ (44) 

n — >-oo 

lim (D m - k f n *(x)D k f n (x))^ 

n — >oo 

where D k f n and D m _ k f n are the k x k and m — k x m — k diagonal blocks of 
the full Jacobian. For k — m , £. = are the Lyapunov exponents. 

Proposed by Pecora & Carroll [1990, 1991] to characterize synchronization in 
chaotic systems, rigorous conditions for the existence of these limits have been 
proven in [Vilela Mendes, 1998]. Existence /i-almost everywhere of both Lya- 
punov and conditional exponents is guaranteed by the conditions of Oseledec's 
multiplicative ergodic theorem, in particular the integrability condition, 

J ii(dx)\og + \\T(x)\\ < oo (45) 

T being either the Jacobian or its k x k and m — k x m — k diagonal blocks. 
The set of points where the limit is defined has full measure and 

lim -log\\D k f n (x)u\\ =c| fe) (46) 

n — >oo fi 

with 7^ u € E l x /E^ +1 , E l x being the subspace of R k spanned by eigenstates 
corresponding to eigenvalues < exp(^| fc - ) ). 
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Based on the spectra of Lyapunov and conditional exponents, several global 
quantities have been defined to characterize self-organization and creation of 
structures in networks of multiagent systems with arbitrary connection struc- 
tures. I list here the definitions and refer to [Vilela Mendes, 2000b, 200f] for 
proofs and examples. 

3.1.1 Structure index related to the Lyapunov spectrum 

A structure (in a collective system) is a phenomenon with a characteristic scale 
very different from the scale of the elementary units in the system. In a multi- 
agent system, a structure in space is a feature at a length scale larger than the 
characteristic size of the agents and a structure in time is a phenomenon with 
a time scale larger than the cycle time of the individual agent dynamics. A 
(temporal) structure index may then be defined by 



where N is the total number of components (agents) in the coupled system, N s 
is the number of structures, Tj is the characteristic time of the structure i and 
T is the cycle time of the isolated agents (or, alternatively the characteristic 
time of the fastest structure) . A similar definition applies for a spatial structure 
index, by replacing characteristic times by characteristic lengths. 

Structures are collective motions of the system. Therefore their character- 
istic times are the characteristic times of the separation dynamics, that is, the 
inverse of the positive Lyapunov exponents. Hence, for the temporal structure 
index, one may write 



the sum being over the positive Lyapunov exponents Aj. Ao is the largest Lya- 
punov exponent of an isolated component or some other reference value. 

The temporal structure index diverges whenever a Lyapunov exponent ap- 
proaches zero from above. Therefore the structure index diverges at the points 
where long time correlations develop. Also, when in a multiagent network the 
coupling between the agents increases, the positive part of the Lyapunov spec- 
trum contracts leading to an effective dimension reduction and to partial syn- 
chronization effects [Vilela Mendes, 1999]. 

3.1.2 Exponent entropies and dynamical selforganization 

Self-organization in a system concerns the dynamical relation of the whole to 
its parts. The conditional Lyapunov exponents, being quantities that separate 
the intrinsic dynamics of each component from the influence of the other parts 




(47) 




(48) 
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in the system, provide a measure of dynamical self organization Is (/z) 

N 

IM = E ^ hk ^ + h ™-k{v) - h(ji)} (49) 
fc=l 

the sum being over all relevant partitions Sfc = R k x R m ~ k and 

w= E d fe) ;^- fe (M)= E el m - fe) ;^)=E A » 

€ f)>o f< m " fe) >o A *>° 

are the exponent entropies, that is, the sums over positive conditional and Lya- 
punov exponents. 

may also be given the following dynamical interpretation: Lyapunov 
exponents measure the rate of information production or, equivalently, they de- 
fine the dynamical freedom of the system, in the sense that they control the 
amount of change that is needed today to have an effect on the future. In this 
sense the larger a Lyapunov exponent is, the freer the system is in that particu- 
lar direction, because a very small change in the present state will induce a large 
change in the future. The conditional exponents have a similar interpretation 
concerning the dynamics as seen from the point of view of each agent and his 
neighborhood [Vilela Mendes, 2000b]. However the actual information produc- 
tion rate is given by the sum of the positive Lyapunov exponents, not by the 
sum of the conditional exponents. Therefore, is a measure of apparent 

dynamical freedom (or apparent rate of information production). 

Being constructed as functions of well defined crgodic limits, both and 
S are also well defined crgodic parameters. They characterize the dynamics of 
multiagent networks and, in addition, also provide some insight on the relation 
between dynamics and the topology of the network [Araiijo et al, 2003]. 



3.2 The problem of pattern discovery: Computational me- 
chanics 

The ultimate practical goal in the study of dynamical systems is the construction 
of models by which these systems might be predicted and (or) controlled to some 
useful purpose. In extended systems with many degrees of freedom, unless exact 
solutions are known, even a knowledge of the (microscopic) equations of motion 
might not be very useful to predict the collective features and patterns that the 
system generates. Crutchfield and collaborators [Crutchficld & Young, 1989] 
[Crutchfield, 1994] [Shalizi & Crutchficld, 2001] have developed a program of 
pattern discovery and construction of minimal models inferred directly from 
the data generated by the dynamical systems. Central to this approach is the 
notion of causal state. Given the knowledge of the infinite past of a system, 
a causal state is an equivalence class of pasts that have the same conditional 
distribution of futures. Denoting by s and s the semi-infinite past and future 
time sequences of coded states of the system, two past sequences si and si 
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belong to the same causal state if 



for all s (except perhaps in a zero measure set). The dynamics of the system 
is then characterized by the set of causal states and the transition probabilities 
between them. That is, the system is mapped into a non-deterministic automa- 
ton called an e— machine. Minimality and uniqueness of the e— machines has 
been proved. Although more general, this approach bears some relation to the 
reconstruction of hidden Markov processes and to grammatical inference. 

As a tool for network dynamics this approach might be useful whenever 
analytical equations are intractable or unknown. Reconstruction algorithms 
for e— machines were developed in some cases [Hanson & Crutchheld, 1997] 
[Crutchfield & Feldman, 1997]. For extended system with many degrees of 
freedom and irregular connections, one problem might be the large size of the 
causal state alphabet. Nevertheless this is a very interesting general approach 
that might be useful to map network dynamics onto probabilistically equivalent 
automata. 

3.3 Construction of invariant measures 

In general a deterministic system has a multitude of invariant measures. How- 
ever, some of them have little practical interest, because they are not stable for 
small random perturbations. Because systems in Nature are subjected to per- 
turbations, only the stable measures are physical measures. In some cases it is 
possible to use the properties of the deterministic system to identify the physical 
measures. For example, in Axiom A systems a unique physical measure may 
be identified with the Sinai-Bowen-Ruelle (SBR) measure, a measure absolutely 
continuous along unstable manifolds. However in most cases, for example in the 
multistable systems so frequent in natural networks, the SBR characterization 
is useless. Instead, it is better to study the stochastic differential equation that 
is obtained from (JIJ by addition of a small noise term 



Wt being a Wiener process and a (X) a A— dependent diffusion coefficient. A 
great deal of information on the invariant measure for this process may be 
obtained using the theory of small random perturbations of dynamical systems 
[Freidlin & Wentzell, 1984, 1994] [Kifer, 1974]. 

If, in the decomposition (JSJ, X (x) has only a gradient component, an explicit 
form for the invariant measure may be obtained. If 



dxi = Xi 0) dt + ea (X) dW t 



(50) 



X{x) 



V (9 )V(x) 



(51) 



V( a ) being the gradient in the metric 




(52) 
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with a (x) in l|5U|) chosen such that 



a-ij (x) = (a (x) a* (x))^ = g tj (x) (53) 

then, the density of the invariant measure is 

P E (x) = C £ exp {-2e~ 2 V (x)) (54) 

as may be easily checked from the forward Kolmogorov equation. In this case, 
finding the stable minima and level sets of V (x) one characterizes the multista- 
bility of the network, their basins of attraction and, from the values of V (x) in 
these sets, the relative occurrence probability of each attractor. 

For general X (x) , small e estimates of the invariant measure for (|50|) are 
also possible. Here the crucial role is played by the functional 



1 f T 

Sot (<p) = 2 J o ai i 



t) <p\ -X* fa) )[<fi-X> (<p t ) dt (55) 



and the infimum 

U (x, y) = inf {S 0T {(f) : ipo = x, ip T = y,t £ [0, T}} (56) 

taken over intervals [0, T] of arbitrary length. 

An equivalence relation is established between points in the domain by x ~ y 
if U (x,y) = U (y,x) = 0. Let the domain be partitioned into a number of 
compacta {Ki] with each oj— limit set of the deterministic dynamics contained 
entirely in one compactum and x ~ y inside each compactum. Then, the (small 
e) asymptotics of the invariant measure is obtained from the invariant measure 
of the Markov chain of transitions between the compacta. For sufficiently small 
e the measure of each compactum is approximated by 

exp j-e~ 2 (w {Ki) - mmW{Ki)\ } (57) 

where 



W{K i )= min \ V{K m ,K n ) (58) 
{m—>n)Eg 

V {K m , K n ) is the minimum of the function (|56|l between points in compacta 
K m and K n and the sum runs over graphs that have exactly one closed cycle 
and this cycle contains the compactum Ki. For proofs I refer to [Freidlin & 
Wentzell, 1984]. 



3.4 A family of ergodic parameters 

Ergodic parameters like the Lyapunov and the conditional exponents, are global 
functions of the invariant measure. However, the invariant measure itself con- 
tains more information. Ergodic parameters being defined by infinite-time lim- 
its, these quantities will fluctuate and, in general, fluctuations will not be Gaus- 
sian. The quantity describing the fluctuations is again an ergodic parameter 
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and the same reasoning applies in turn to its fluctuations, etc. [Ruelle, 1987]. 
Therefore, to characterize the measure, a larger set of parameters is needed. 
To construct this larger set from the fluctuations is not very practical and a 
different approach will be followed here, namely a variational approach. 

In a restricted sense, a variational principle states that the equations of 
motion may be written in the form SS = 0, where S is a functional of the 
dynamical variables and S is the Gateaux derivative. Only a limited set of 
dynamical systems may be described by a variational principle in this restricted 
sense. However, if one only requires that SS = and the equations of motion 
possess the same set of solutions, essentially all differential equation problems 
admit a variational formulation [Tonti, 1992]. Let 

Xi= Xi (x) (59) 
be a differentiable continuous-time dynamical system and S be the functional 



S = j Jo dtdT Yl{^ (*) (*(*)) }ff(*> T ){^ (t)-^O(t))} 



(60) 



where g(t,r) is a symmetric kernel (g(t,r) — g(r,t)). Let us compute the 
Gateaux derivative for variations in space restricted by the boundary conditions 

u (0) = u (T) = (61) 

From 

S U S = -j£ dtdr «*(*) {<^^f^ - d ^ X i ( x (*)) 9 (t, r)| { k (r) - X, (x 

(62) 

we have 

Lemma: The equations of motion (|59|l and the critical points of the func- 
tional (S U S — 0) have the same set of solutions if 

K (t, t) = hi ^^ ~ d k Xi (x (t)) g (t, t) (63) 

is invertible. 
Remarks: 

a) If K (t, t) is not invertible, the solutions of the equations of motion are 
still critical points of the functional, but this one might have other solutions. 

b) A variational principle, with only u (0) = being required, may also be 
obtained by choosing a kernel such that g (t, T) = 0. 

The critical points of the S functional contain the same information as the 
equations of motion. Therefore the dynamics may be characterized by the prop- 
erties of the critical points, in particular by their Hessian matrix. Computing 
the second Gateaux derivative on the orbits one obtains 

S lvS\ ss=0 =11 dtdr £ m (t) Vj (r) Ha (t, r) (64) 



26 



with 

(65) 

Now assume that the symmetric kernel g (t, r) is a function of finite support of 

t - T 

g (t, T )=g(t-T)=0 for \t - r| > r (66) 



Define 



as well as 



t(™)' _ j(n) 
J 0,T ~ J 0,T 



Then 



jW = y jf Tr (if" (i, t)) <Mt 



(67) 



f + |Tr( J ff n (t,T))|dtdr+ / / |Tr(if n (i,r))|*dr 

T-r J J-r 



t(™)' <- t(«)' i t(™)' 
" , 0,Ti+T 2 — J 0,Ti "T" J T!,T 2 



and we are in the conditions of Kingman's sub-additive ergodic theorem. Taking 
limits, if both Xi and diX^ are bounded J^'and Jq^ differ only by a finite 
quantity and one concludes: 

Theorem: If /i is an invariant measure of the dynamics in (|59l) . Xi and diXk 

(n) 

are bounded and there is M > such that Jg ^ > — M for sufficiently large T, 



then the limit 



exists and 



I n ( M ) = lim i jg (68) 



J — >oo 1 ' J — >oo 1 I ' 



/„ (//) for n = 1, 2, • • • is a family of ergodic parameters for the fi— measure 
preserving dynamics. 

A similar construction for discrete-time maps may be found in [Vilela Mendes, 
1984] [Carreira et at, 1991]. 

3.5 Synchronization, mode- locking and dynamical corre- 
lations 

The onset of correlated motions in coupled many-agent systems is a phenomenon 
of widespread occurrence in many scientific fields. The most dramatic effect is 
the synchronization of assemblies of coupled dynamical systems which, when in 
isolation, may have quite different rhythms [Pikovski et at, 2001]. Examples 
are biological rhythms [Winfree, 1967] like the pacemaker cells in the heart [Pe- 
skin, 1975], neural systems [Golomb & Hansel, 2000], synchronous metabolism 
[Aldridge & Pye, 1976], flashing fireflies [Buck, 1988], laser arrays [Jiang & Mc- 
Call, 1993], even fads and social trends may be interpreted as synchronization 
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of distinct agent dynamics. The study of the correlated behavior of many-agent 
dynamics is also closely related to the problem of control in extended dynamical 
systems. 

I will consider both the coupled behavior of non-chaotic systems (oscillators 
with distinct individual frequencies) and of systems with isolated chaotic dy- 
namics. In both cases one may distinguish between globally coupled systems 
and systems where each agent has a limited range or number of interacting 
partners. 



Kuramoto (K=y) 




20 40 60 
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Figure 11: A Kuramoto system below threshold 

For systems of oscillators the canonical example is the Kuramoto model 
[Kuramoto, 1984, 1991], 



dt 



K 



N - 



N 



sin ( 



(69) 



with K > and the frequencies d>i randomly distributed around a central value 
loo with the shifted Cauchy distribution 



7 



7 2 + (uj - uj y 



(70) 



A great deal of work has been done on this model (for a review see Strogatz 
[2000]). The existence of a synchronized cluster is characterized by the order 



28 



Kuramoto (K=5y) 
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Figure 12: A Kuramoto system above threshold 



parameter 



•(*) = 



1 N 

-Y 



(71) 



It is found that in the N — > oo and t — > oo limit, r — for K < 27 and 
r = v/1 — (2'y/K) for X > 27. That is, there is a coupling threshold above 
which part of the oscillators starts to synchronize. Figs. 11 and 12 show the 
nonsynchronizcd (at K = 7) and the synchronized (at K — 57) behavior for 
100 oscillators. The upper plot displays the color-coded values of the oscillator 
variables at the end of each unit time interval. The lower plots show the time 
evolution of the order parameter. Fig. 13 compares the numerically computed 
Lyapunov spectrum in the synchronized and non-synchronized situations. One 
sees that even below the synchronization threshold (K = 2j) , part of the Lya- 
punov exponents becomes negative, meaning that there are many contracting 
directions, implying an effective dimension-reduction of the asymptotic behavior 
of the system. This clearly suggests that synchronization is not the whole story 
and that even before synchronization strong correlations must develop between 
the dynamics of the individual oscillators. 

A type of correlation, of which synchronization is a limiting case is mode- 
locking. Mode-locking is the entrainmcnt of some integer combination of the 
frequencies to zero. It also plays an important role in the dynamics of cou- 
pled oscillators [MacKay, 1994]. However even if all the effective frequencies are 
incommensurable, the existence of negative Lyapunov directions, implies the 
existence of dynamical correlations between the oscillators. What is important 
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Figure 13: Lyapunov spectrum below and above threshold for the Kuramoto 

system 

is the dimension of the invariant measure and the correlations may be charac- 
terized by the eigenvectors of the Lyapunov spectrum. These notions are better 
clarified in a simple model with exactly computable Lyapunov spectrum. Let 
the dynamics of an assembly of discrete-time oscillators be 

k A 

Xi (t + 1) = Xi (t) + Ui + — — - 22 fa {xj - %i) (72) 

i=i 

with Xi € [0, 1) and f a (xj — Xi) = a (xj — Xi) (modi) and the lo[s distributed 
according to p (u), as above. 

The Lyapunov spectrum is composed of one isolated zero and log (l - ^ak 

(N — l)-times. However, although (for all k > 0) N — 1 contracting directions 
are always present, it is only for sufficiently large k that synchronization effects 
emerge as shown in Figs. 14 and 15. Nevertheless dynamical correlations do exist 
for all k, no matter how small and the Lyapunov dimension is always one. In 
this case, the eigenvectors of the Lyapunov spectrum may be exactly computed 
and the correlations explicitly identified. This is illustrated in Fig. 16. 

So far I have dealt with coupled oscillators, that is, with systems which have 
individual nonchaotic dynamics. Another important field with many practical 
applications refers to the case where the individual node dynamics is chaotic. 
Synchronization of chaotic systems has been extensively studied (for a review 
see Pecora et al. [1999]) and is still a field of current research [Wei et ai, 
2002]. However, as in the oscillators, for networks of chaotic elements the in- 
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Figure 14: Non-synchronized behavior of the discrete-time oscillators (Ea l72|l 

teresting phenomena go beyond synchronization effects. Correlations and self- 
organization effects may be characterized by ergodic parameters. I refer to 
[VilelaMendes, 1999, 2000b] for illustrative examples of networks of chaotic el- 
ements both globally connected and with a limited range of interactions. The 
Lyapunov spectrum and the entropies associated to the conditional exponents 
provide a characterization of the emergent phenomena. It should be noticed 
that dynamical correlations play an important role in the organization of the 
dynamics even when there is no reduction of the Lyapunov dimension [Vilela 
Mendes, 1999]. As before these correlations are associated to the eigenvectors 
of the Lyapunov spectrum. 

3.6 Dynamics and network topology 

The topology of the network connections has a determining effects on the dy- 
namical phenomena taking place in the network [Watts, 1999]. Local clustering 
and far-reaching connections in the network influence the spread of infectious 
diseases [Boots & Sasaki, 1999] [Keeling, 1999] [Pastor-Satorras & Vespignani, 
2002] [Lloyd & May, 2001] or social fads. The nature and range of the couplings 
influences the speed of signal propagation and even the computational abilities 
of the network. On the other hand the topological structure of the networks by 
influencing the dynamics, may have a feedback effect on network growing and 
therefore also evolutionary significance. 

In particular the small world topology [Watts & Strogatz, 1998] [Watts, 
1999] [Strogatz, 2001] (with both small path length and large clustering) has 



31 



Syncnet (k=0.8 




'~Y~ J — — yv — \[^p{T^fYYTTT^ — — — ' — ' — — V~Y 



100 120 140 160 



Figure 15: Synchronized behavior of the discrete-time oscillators (Ea 172)1 



been found to modify or enhance coherent behavior effects [Lago-Fernandez 
et al., 2000] [Gade & Hu, 2000] and in general influence the dynamics in the 
network [Yang, 2001] [Kulkarni et al., 1999]. An attempt has also been made 
to relate the ergodic parameters of a dynamical system, living on the network, 
with the changes of topological structure associated to the small- world features. 
It turns out that whereas the emergence of short path lengths is associated to a 
modification of the Lyapunov spectrum, the transition from the small world to 
the random regimen is characterized by the conditional exponent entropies (for 
details see Araujo et al. [2003]). 



4 The logic approach to network dynamics 

Thomas and collaborators [Thomas & D'Ari, 1990] [Thomas, 1991] [Thieffry et 
al., 1993] [Snoussi & Thomas, 1993] [Thomas et al., 1995] have developed logical 
tools to analyze biological regulatory networks. This treatment seems particu- 
larly appropriate to deal with regulatory networks where most interactions are 
characterized by a threshold and a saturation plateau. The regulator is usually 
inefficient below a threshold concentration and its effect rapidly levels off above 
the threshold. 

The elements of the network and their interactions are represented by dis- 
crete variables, functions and parameters. Because some variables have several 
distinct actions, one often needs to consider more than two logical levels. One 
may also consider the threshold values s^' separating the logical values. Then, 
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Figure 16: Correlations in the discrete-time oscillators system. 

the most general logical variable takes values in the set 

{0, S « 1, S ( 2 \2,...} (73) 

The state of the system is described by the state vector xyz ■ ■ ■ containing the 
logical values of the variables and the evolution of the system by a vector of 
functions XYZ ■ ■ ■ representing the excitatory or inhibitory interactions in the 
network. For example, for the graph of interactions 

x — ► y — ► z (74) 



the logical functions are 

X=z, Y = x, Z =y (75) 

X = means that the product x is not being produced and X = 1 means that 
the product x is being produced, not that x = 1 immediately. For example, if 
the initial state is 000, after a certain time the state might become 100 or 001, 
depending on whether the time delay (t x ) for production of x is smaller or larger 
than the time delay (t z ) for the production of z. Notice that although the image 
XY Z of the state 000 is 101, the next state is either 100 or 001 because it is 
highly unlikely that t x = t z . The state 001 is a stable state because it equals 
its image, whereas 100 further evolves to 110 or 101, of which 110 is stable but 
101 is not, etc. 
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A central role in this formalism is played by the oriented circuits (closed series 
of interactions where each element appears only once), called feedback loops. 
Feedback loops are positive or negative according to whether they have an even 
or odd number of negative interactions. Positive loops generate multistability 
and negative loops generate homeostasis, that is, the variables in the loop tend 
to middle range values, with or without oscillations. 

The role of the parameters in the logical approach is to allow for the dis- 
tinction between weak and strong interactions. Therefore the parameters are 
actually real values, not logical variables. For example, in the following inter- 
action graph 

x +±± y □+ (76) 

x has two values (0, 1) because it has one action only and y has three values 
(0, 1, 2) because it has two actions and thus two thresholds. The corresponding 
logical functions are 

X=y\ Y = x 1 +y 2 (77) 

in which all the variables are Boolean. The function Y means that if x is above 
its threshold or y above the second threshold (x 1 = 1 or y 2 = 1) then y is going 
to be produced. However we might give different weights to the variables by 
writing 

X = d x ^Ki y 1 ^, Y = d y (Kzx 1 + K 3 y 2 ) (78) 

where the K i s are real variables and d x and d y operators that discrctizc the 
value in brackets. Here the + sign is the real algebraic sum not the logical one. 

States with variable values xyz ■ ■ ■ that involve only logical values are called 
regular states and those which involve some threshold values are called singular 
states. It is found that each feedback loop can be characterized by a singular 
logical state located on the thresholds. For appropriate parameter values, this 
state is stationary and the corresponding loop functional. In this context, func- 
tional means that if the loop is positive it actually produces multistability and 
if negative it generates homeostasis. A technique for the analysis of the network 
consists in dissociating it into its feedback loops and checking the dynamics of 
each loop. 

Other logical approaches have been developed to analyze networks: the 
Boolean network models [Kaufman, 1993] [Somogyi & Sniegoski, 1996], hy- 
brid models using logical and continuous variables [Glass, 1975] [Lewis & Glass, 
1991] and rule-based formalisms [Brutlag et at, 1991] [Meyers & Friedland, 
1984] [Shimada et al, 1995] [Hofestadt & Meineke, 1995]. 

A general problem in the logical approaches to network dynamics is to es- 
tablish the correspondence of the logical description to the corresponding set of 
differential equations. For the multilevel approach (with parameters) of Thomas 
and collaborators, if the thresholds are represented by Hill functions, the cor- 
respondence becomes perfect when the Hill functions become very steep. The 
problem of establishing the correspondence between logical descriptions and 
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smooth dynamical systems has also been addressed in other contexts. For ex- 
ample, in [Martins & Vilela Mendes, 2001] a correspondence is established be- 
tween a type of neural networks and a logical programming language and in 
[Martins et at, 2001] a correspondence between controlled dynamical systems 
and context-dependent languages. Similar techniques may be used for general 
networks [Aguirre et al., 2004]. 
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